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The  progressive  elevation  of  clock  speeds  in  digital  electronic  circuits  is  compelling  the  merging  of  design  considerations  for 
such  circuits  with  those  of  analog  microwave  circuits.  For  sufficiently  high  clock  rates  digital  circuits  develop  substantial 
analog  effects,  such  as  coupling  and  radiation,  that  markedly  affect  the  flow  of  signal  pulses  in  such  structures  as  multilayer 
circuit  boards  and  multichip  modules  resulting  in  unreliable  circuit  or  device  operation.  Traditional  circuit  models  such  as 
SPICE  cannot  account  for  the  distributed  field  physics  inherent  in  these  processes.  To  accurately  analyze  and  design 
complex  digital  circuits  above  500  MHz,  it  is  necessary  to  resort  to  full-wave  Maxwell’s  equations  simulations  to  rigorously 
compute  the  electromagnetic  field  environment  wifliin  the  circuit.  The  goal  of  this  effort  was  to  initiate  the  development  of 
electromagnetic  models  which  serve  as  a  bridge  between  electromagnetic  field  analysis,  nonlinear  circuit  analysis  and 
CAD/CAM  design  tools  in  order  to  model  high  speed  digital  circuits.  The  computational  methodology  is  based  on  a  finite 
difference-time  domain  (FDTD)  solution  of  Maxwell’s  equations  which  lends  itself  to  analyzing  digital  circuits  and  devices 
accurately  in  the  presence  of  complex  geometry  and  materials  and  general  digital  signals.  The  bridging  models  permit  the 
results  of  FDTD  simulation  to  be  formulated  in  terms  of  circuit  q[uantities  that  can  be  used  as  input  into  circuit  programs  such 
as  SPICE.  Such  circuit  quantities  now  contain  the  mutual  coupling  effects  inherent  in  the  full  wave  Maxwell’s  equation 
solution  making  the  circuit  solution  much  more  accurate. 
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Electromagnetic  Field  Analysis  of 
Digital  Circuits  and  Modules  for 
Electromagnetic  Compatibility  and  Reliability 


1.0  OBJECTIVE  AND  SCOPE 

The  goal  of  this  project  was  to  initiate  the  development  of  bridging  software  between 
electromagnetic  field  analysis,  nonlinear  circuit  analysis,  and  CAD/CAM  tools  to  model  high¬ 
speed  digital  circuits  and  multichip  modules  for  electromagnetic  compatibility  and  reliability, 
"^e  computational  methodology  was  based  upon  the  finite-difference  time-domain  (FD-TD) 
solution  of  Maxwell's  equations.  The  result  of  this  research  should  be  the  initiation  of  a  novel 
computational  simulation  tool  that  would  assist  understanding  the  operation  of  high-speed 
electronic  devices  of  high  potential  value  to  Rome  Laboratory. 

This  was  initially  a  six-month  program,  later  extended  at  no  additional  cost  to  eight  months, 
that  provided  partial  support  for  one  Ph.D.  graduate  at  Northwestern.  Assistance  in  running 
FD-TD/SPICE  simulations  was  provided  at  no  cost  by  Prof.  Melinda  Piket  May  and  her  students 
at  the  University  of  Colorado  at  Boulder. 

2.0  BACKGROUND 

The  progressive  elevation  of  clock  speeds  in  digital  electronic  circuits  is  compelling  the  merging 
of  design  considerations  for  such  circuits  with  those  of  analog  microwave  circuits.  This  is 
because  digital  circuits  develop  substantial  analog  effects  when  their  clock  rates  are  high  enough. 
These  effects  include  such  fuU-wave  electromagnetic  field  phenomena  as  coupling  and  radiation, 
and  can  markedly  affect  the  flow  of  signal  pulses  in  typicfd  structures  such  as  multilayer  circuit 
boards  and  multichip  modules. 

In  fact,  coupling  and  radiation  can  cause  unreliable  circuit  operation.  Traditional  circuit 
models  such  as  SPICE  do  not  have  the  capability  to  calculate  the  distributed-field  physics 
inherent  in  these  processes.  Especially  for  digital  circuits  operating  above  500  MHz,  it  may  not 
be  possible  to  insure  reliable  operation  without  resorting  to  full-wave  Maxwell's  equations 
simulations.  Such  simulation  tools  must  properly  account  for  the  physics  of  UHF  and 
microwave  electromagnetic  wave  energy  transport  along  metal  surfaces  such  as  ground  planes, 
or  in  possibly  inhomogeneous  dielectric  media  away  from  the  metal  signal  paths.  Such  tools 
must  also  be  capable  of  dealing  with  the  complexity  of  modem  engineering  designs. 

A  strong  candidate  for  this  role  is  FD-TD,  which  is  currently  in  worldwide  use  in 
applications  relating  to  radiated  and  scattered  electromagnetic  wave  interactions  with  material 
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structures.  This  report  discusses  the  application  of  FD-TD  modeling  to  analyze  the 
electromagnetic  wave  properties  of  high-speed  electronic  circuits  in  pneral,  making  no 
distinction  between  digital  or  analog.  We  consider  the  analysis  of  both  passive  interconnects  and 
circuit  devices,  including  extraction  of  lumped-circuit  equivalences,  5-parameters,  and 
connections  to  linear  and  nonlinear  loads. 

3.0  TASK  1  ■  EXTRACT  LUMPED  DEVICE  MQPELS 
FROM  Fn-TD  ELECTROMAGNETIC  FIELD  DATA 

This  section  reports  means  to  extract  equivalent  transmission  line,  inductance,  and  capacitance 
parameters  from  FD-TD  simulations  of  general  3-D  electrical  structures.  The  structures  include 
interconnects,  vias,  connectors,  and  power-distribution  systems. 

The  analysis  begins  with  properly  connecting  the  computed  electric  and  magnetic  field 
vector  distributions  in  time  and  space  to  the  usual  circuit  quantities  of  voltage  and  current: 

V(/,x,.)  =  f  I(uxi)  =  ^H(ux)‘dl  (1) 

Here,  Q  is  a  contour  extending  from  a  defined  voltage  reference  point  (usually  a  ground  plane) 
to  the  circuit  at  location,  x,..  In  many  cases,  Xj  is  a  point  on  a  metallic  strip  transmission  line  (a 
microstrip)  that  propagates  the  dominant  TEM  mode.  In  this  situation,  V(r,x,  )  is  independent  of 
the  choice  of  Q  if  this  path  is  confined  to  the  transverse  plane,  and  can  be  conveniently 
chosen  to  extend  in  a  perpendicular  manner  from  the  ground  plane  to  the  adjacent  surface  of  the 
microstrip.  Similarly,  the  selection  of  the  contour  Q  to  wrap  completely  around  the  strip 
conductor  at  its  surface  in  the  transverse  plane  provides  the  local  current.  Care  must  be  exercised 
in  applying  (1)  if  there  is  evidence  that  the  circuit  path  of  interest  is  propagating  non-TEM 
modes;  the  accuracy  of  this  formulation  can  degrade  rapidly  if  the  TEM  assumption  is  not  met. 

3.1  Transmission  Line  Parameters  and  Load  Impedance 

Wideband  frequency-domain  transmission  line  parameters  can  be  found  by  applying  the  Fourier 
transform  to  the  voltage  and  current  responses  of  (1)  for  an  impulsive  excitation  of  the  line  [1]. 
For  example,  the  line  impedance  as  a  function  of  frequency  is  obtained  from 

Zo(m,Xi)  =  :r[v(f,x,.)]/^/(Ax..)]  (2) 

where  ir[  ]  is  the  Fourier  transform  operator.  Further,  if  gC^.x.)  denotes  a  voltage  or  current 
waveform  at  x  =  x,-,  and  g(^xp  is  the  corresponding  waveform  at  x  =  x^,  then  the  propagation 
constant  y  can  be  calculated  using  the  following  relationship: 

e-'W-  (3a) 

where  d  =  x^  -Xj.  Rearranging  this  expression  gives  y  as  a  function  of  frequency: 
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(3b) 


7(0)  =  ^In-I 


!F 

9 

Defining  7(0)  =  a(o)+ jp((o),  the  group  velocity  v^(o)  is  calculated  as 


■.«  -  [^] 


rl 


(4) 


For  an  impulsive  line  excitation,  the  wideband  reflection  coefficient  from  a  given 
transmission  line  load  is  calculated  using 


r(Q>,x,.) 


7 

7 

V,(l.x,) 

(5a) 


where  V,  is  the  reflected  voltage  and  is  the  incident  voltage  observed  at  x,..  This  reflection 
coefficient  is  subsequently  transformed  to  the  plane  of  the  load  via 

ri(o)  =  no.x.)  (5b) 


where  £  is  the  distance  from  x,.  to  the  load.  The  effective  load  impedance  is  then 


Zjo)  = 


Zo 


i+r,(o) 

i-r,(o). 


(6) 


where  Z^(o),  ;^(o,x,),  and  r^(G))  are  all  complex  values. 

The  discussed  derivation  of  impedance  does  not  take  into  account  that  the  voltage  and 
current  values  derived  from  corresponding  electric  and  magnetic  fields  in  the  Yee  grid  are  offset 
from  each  other  by  one-half  space  cell  and  one-half  time  step.  Ignoring  these  offsets  can  lead  to 
errors  for  certain  geometries.  Reference  [2]  reported  a  simple  but  effective  interpolation  that 
permits  the  voltage  and  current  data  used  in  the  impedance  calculation  to  be  at  the  same  space- 
time  point  Given  a  voltage  Vj.  ^  *  and  adjacent  currents  4+i/2,7.t  desired 

spatially  interpolated  current  is 


'4 


“ 


4-1 


li+l/2,7.t  'li-1/2,;.* 


(7) 


The  half  time  step  can  be  accounted  for  by  multiplying  Vj,.  *  by  a  factor  of 
corrected  equation  for  impedance  now  becomes 


(8) 
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Reference  [3]  reported  studies  to  validate  FD-TD  modeling  of  line  impedance  for  the 
important  case  of  parallel  coplanar  microstrips.  The  first  studies  were  done  by  setting  up  FD-TD 
grids  for  single  x-directed  microstiips  of  negligible  metalization  thickness  and  a  variety  of  widths 
over  dielectric  substrates  on  the  order  of  1  mil  (0.001  inch)  thick.  Grid  cell  size  A  for  these 
cases  was  on  the  order  of  0. 1  roil  with  At  on  the  order  of  4.2  fs.  In  all  cases,  conductors  were  on 
grid";  that  is,  they  were  located  at  planar  loci  of  tangential  E  components  in  the  FD-TD  lattice 
that  were  set  to  zero  for  all  time  steps.  Excitation  of  a  microstrip  was  provided  by  specifying  a 
Gaussian-pulse  time  history  for  a  group  of  collinear  electric  field  components  (usually  E^) 
bridging  the  gap  between  the  ground  plane  and  the  strip  conductor  at  the  desired  source  location. 
Application  of  (2)  showed  the  FD-TD-computed  values  of  Zq  to  be  virtually  independent  of 
frequency  up  to  1.0  GHz,  and  on  the  order  of  1%  agreement  with  textbook  values  [4]. 

Next,  using  similar  FD-TD  grid  resolutions,  [3]  considered  single  microstrips  with  finite 
metalization  thickness,  possibly  fully  embedded  within  a  dielectric  layer.  Here  FD-TD 
predictions  were  compared  to  measurements.  In  one  example,  a  1.1-mil-thick,  1.4-mil-wide 
metal  strip  was  assumed  to  be  suspended  1.1  mils  above  a  large  ground  plane  within  a  3.1 -mil- 
thick  dielectric  layer  having  £,.=3.2.  The  FD-TD  numerical  simulation  predicted  a  flat 
characteristic  impedance  of  48  up  to  1.0  GHz,  agreeing  with  the  experimental  results  to  within 
the  measurement  uncertainty  (about  0.2  £2).  The  computed  variation  of  Zq  above  1  GHz  was 
found  to  be  ±  2  £2.  Similar  excellent  agreement  was  found  for  the  propagation  delay,  where  the 
experimental  value  was  150.5  ±  1.5  ps/inch,  while  the  FD-TD  prediction  was  149.5  ps/inch. 

Reference  [3]  also  reported  a  study  of  the  impedance  and  propagation  delay  of  three  parallel, 
coplanar  microstrips  having  finite  metalization  thickness.  Each  microstrip  had  the  geometry 
discussed  above  and  was  separated  by  3.6  mils  from  the  adjacent  lines.  Even-mode  results  were 
obtained  with  all  three  strips  excited  simultaneously  with  the  same  polarity,  while  odd-mode 
results  were  obtained  with  the  two  outer  strips  excited  with  the  opposite  polarity  relative  to  the 
center  strip.  Results  were  also  observed  for  the  center  strip  excited  with  the  two  outer  strips 
floating.  The  FD-TD  modeling  showed  that  signal  propagation  on  adjacent  lines  significantly 
affected  the  impedance  of  the  line  of  interest,  and  to  a  lesser  degree  affected  its  propagation 
delay.  These  results  were  confirmed  by  laboratory  studies  which  showed  an  approximate  7  £2 
elevation  of  the  characteristic  impedance  for  the  even-mode  excitation,  and  an  approximate  7  £2 
reduction  of  the  impedance  for  the  odd-mode  excitation  (both  from  dc  to  about  1  GHz). 

3.2  .y-Parameters 

Given  a  multiport  network,  wideband  complex-valued  scattering  parameters  S^„  can  be  obtained 
for  an  impulsive  excitation  as  follows  [1,5]: 


(CD  X  X  )  = 


(9) 


where  V„(o),x„)  is  the  voltage  at  port  m  at  observation  plane  x„,  V„(,o),x„)  is  the  voltage  at  port 
n  at  observation  plane  x„,  and  and  Zo  „  are  the  characteristic  impedances  of  the  lines 
connected  to  these  ports.  For  example,  to  obtain  5i,,  the  incident  and  reflected  pulses  at  port  1 
must  be  known.  To  obtain  S^i,  we  must  observe  the  transmitted  pulse  emerging  at  port  2 
corresponding  to  the  known  incident  pulse  at  port  1. 


In  most  cases,  the  magnitudes  of  the  5-parameters  are  the  primary  data  used  by  engineers  to 
characterize,  for  example,  the  filtering  properties  of  a  network.  The  magnitude  data  are 
independent  of  the  observation  positions  on  the  transmission  lines  feeding  the  corresponding 
ports,  assuming  that  the  feeding  lines  are  either  infinitely  long  or  matched  at  their  far  ends. 
However,  the  phases  of  the  5-parameters  are  clearly  a  function  of  the  positions  of  the  observation 
planes.  As  discussed  later,  phase  data  can  be  important  when  extracting  a  lumped-circuit 
equivalent  network  from  the  observed  5-parameter  variation  with  frequency  [6]. 


3.3  Differential  Capacitance  Calculation 

The  capacitance  per  unit  length  of  a  strip  line  can  be  calculated  indirectly  using  (2)  and  (4)  as 


Cid)  = 


1 

v,(n>)Zo(Q)) 


(10) 


Of  more  interest  is  the  possibility  that  this  capacitance  can  also  be  calculated  directly  from  the 
electric  field  that  FD-TD  simulations  provide.  Given  a  designated  incremental  Ax  Ay  surface 
patch  AA  on  the  bottom  of  a  circuit  trace  (facing  the  ground  plane).  Gauss's  Law  can  be  used  to 
calculate  the  electric  charge  A0  on  this  patch  from  the  normal  component  of  E  originating 
from  the  patch: 


Ae  =  (lla) 

Here  D  =  eE  and  5^  is  a  virtual  surface  positioned  parallel  to  the  trace  and  between  AA  and  the 
ground  plane  that  captures  all  of  the  electric  flux  emanating  from  AA .  Combined  with  the 
voltage  expression  of  (1),  the  incremental  trace  capacitance  (in  farads/meter)  is  given  by  the  ratio 
of  the  charge  on  AA  to  the  voltage 


Ax  V 


(lib) 


This  approaches  the  true  differential  capacitance  per  unit  length  in  the  limit  that  Ax  «  A  for  the 
smallest  relevant  wavelength  A  present  on  the  trace.  In  this  quasi-static  limit,  AC/ Ax  is 
independent  of  the  precise  nature  of  AA ,  5^,  and  the  excitation  of  the  trace. 

Note  that  for  building  physically  representative  models,  this  procedure  lends  itself  to 
calculating  the  capacitance  between  any  two  surface  regions  of  any  shape.  Patch  AA  could  be 
deHned  to  wrap  around  the  entire  circuit  trace.  The  conductor  at  the  other  end  of  the  integral 
path  used  to  deHne  the  voltage  (i.e.,  the  other  node  of  the  capacitor),  could  be  an  adjacent  trace 
rather  than  the  ground  plane.  This  electric  flux  collection  method  can  directly  provide  the 
equivalent  mutual  capacitance  for  subsequent  SPICE  calculations  of  crosstalk. 

3.4  Differential  Inductance  Calculation 

In  a  manner  analogous  to  that  discussed  above,  the  differential  inductance  of  a  circuit  trace 
geometry  can  be  calculated  directly  from  the  magnetic  field  vector  components  that  FD-TD 
simulations  provide.  It  is  well  known  that  the  self  and/or  mutual  inductance  of  any  physical 
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structure  is  entirely  determined  by  the  geometric  relationship  between  the  electric  circuit  that 
transports  the  current  I  and  the  surface  S„  through  which  the  magnetic  flux  generated  by  I 
penetrates.  In  fact,  inductance  is  defined  to  be  the  ratio  of  the  magnetic  flux  penetrating  S„  to  its 
generating  current 


L 


(12) 


To  obtain  inductance  from  the  full-wave  electromagnetic  field,  it  is  crucial  to  properly 
define  the  magnetic  flux  collection  virtual  surface,  S^.  Consider  the  case  of  no  discontinuities  of 
Circuit  Trace  No.  1,  which  carries  the  signal  current,  /„  and  the  underlying  ground  plane,  which 
carries  the  return  current,  Here,  we  define  a  planar  rectangular  "window  pane"  surface,  S„  ^y 
that  has  a  perimeter  bounded  on  the  top  by  Trace  No.  1,  bounded  on  the  bottom  by  the  ground 
plane,  and  extends  a  short  distance  Ax  ^ong  the  trace.  The  incremental  self-inductance  of  Trace 
No.  1  (in  henrys/raeter)  is  given  by 

AAi  ^  Ml  = 

Ax  /i  /, 


where  ^  is  the  magnetic  field  generated  by  /,.  The  incremental  mutual  inductance  relative  to  a 
possible  adjacent  trace  (Trace  No.  2)  is  found  by  obtaining  the  magnetic  flux  penetrating  y  due 
to  current  /j  flowing  on  Trace  2: 

Msl  =  ^21  =  (14) 

Ax  /j  li 

where  Bi  is  the  magnetic  flux  generated  by  /j.  The  inductance  values  of  (13)  and  (14)  approach 
the  true  differential  inductance  per  unit  length  in  the  limit  that  Ax  «X  for  the  smallest  relevant 
wavelength  A  present  on  the  trace.  In  this  quasi-static  limit,  AL/  Ax  is  independent  of  the 
precise  nature  of  and  the  excitation  of  the  trace. 


3.5  Lumped  Inductance  Due  to  a  Discontinuitv 

A  key  element  of  current  engineering  practice  in  the  design  of  high-speed  circuits  is  the 
development  of  lumped-circuit  equivalences  (especially  inductances)  for  discontinuities  in  the 
signal  and  ground  return  paths.  The  equivalent  circuits  for  the  discontinuities  are  then 
substituted  into  SPICE  or  SPICE-like  circuit  modeling  software  to  obtain  the  desired  overall 
circuit  response.  Using  the  Fourier-transformed  FD-TD  fields,  the  desired  equivalent  circuits 
may  be  derived  in  two  ways  [6]:  (1)  the  basic  flux/current  definition  of  inductance,  and  (2)  the 
fitting  of  an  equivalent  circuit  to  the  calculated  variation  of  impedance  or  .S-parameters  with 
frequency.  These  approaches  are  now  summarized. 
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3.5.1  Flnx/Current  Definition 

The  simple  magnetic  flux  collection  procedure  discussed  in  Section  3.4  may  be  extended  to  any 
three-dimensional  circuit  trace  geometry,  including  the  interesting  cases  of  a  via  present  in  the 
signal  path  and  a  geometrical  discontinuity  present  in  the  current  return  path.  Consider  for 
example  a  microstrip  that  traverses  an  open  hole  or  slot  discontinuity  in  the  ground  plane.  Here 
the  use  of  the  planar  "window-pane"  integration  surface  ,  between  the  slot  and  the  circuit 
trace  would  underreport  the  inductance  of  the  discontinuity,  because  a  portion  of  the  magnetic 
flux  in  the  open  hole  or  slot  would  slip  under  its  collection  area.  But  we  note  that  the  magnetic 
flux  is  always  contained  between  the  signal  and  the  return  current  paths.  Thus,  it  is  essential  that 
the  edges  of  S„  are  bounded  by  the  signal  and  return  currents  to  collect  all  of  the  magnetic  flux. 
Then  the  total  magnetic  flux  penetrating  would  be  independent  of  the  precise  surface 
definition  of  5„. 

We  see  that  the  problem  reduces  to  finding  the  return  current  path  and  then  constructing  an 
integration  surface  bounded  by  this  path  and  by  the  signal  trace.  One  possibility  for  is  a  box¬ 
like  or  "hockey  net"  integration  surface.  Similar  to  the  "window  pane,"  this  S„  has  a  perimeter 
bounded  on  the  top  by  the  circuit  trace  of  interest.  However,  this  is  bounded  on  the  bottom 
by  the  ground  plane  along  a  path  that  is  curved  around  the  hole  or  slot  to  follow  the  distorted 
return  current  path.  In  this  manner,  captures  substantially  all  of  the  local  magnetic  flux  and 
properly  yields  the  inductance  of  the  discontinuity.  For  simplicity,  we  assign  four  flat  rectangular 
walls  to  define  the  "hockey  net,"  noting  again  that  the  flux  calculation  is  dependent  upon  the 
location  of  the  edges  of  the  net  rather  than  its  surfaces.  Of  course,  this  Cartesian  integration 
surface  is  natural  for  the  Yee  grid. 

3.5.2  Fitting  270))  or  Sica)  to  an  Equivalent  Circuit 

A  second  approach  for  obtaining  the  equivalent  lumped  inductance  due  to  a  discontinuity  such  as 
a  via  or  a  slotted  ground  plane  is  to  fit  the  calculated  impedance  or  iS-parameter  variation  with 
frequency  to  that  of  an  equivalent  circuit  over  a  frequency  range  of  interest  [7-9].  For  relatively 
complex  geometries  consisting  of  several  package  planes,  vias,  and  microstrips  embedded  in 
multilayered  inhomogeneous  media,  the  fitting  process  can  be  automated  to  an  extent  by  using  a 
software  package  such  as  EEsofs  TOUCHSTONE®  to  manipulate  the  components  in  an 
equivalent  circuit  composed  of  a  number  of  inductor-capacitor-resistor  (LCR)  sections  [9].  The 
inductors  model  the  vias,  microstrip  traces,  and  ground  return  discontinuities,  while  the 
capacitors  account  for  the  coupling  between  the  package  planes.  Resistors  are  associated  with 
internal  impedances  of  sources,  lossy  materials,  and  leakage  radiation  from  the  circuit  planes. 
The  distributed  nature  of  these  processes  is  modeled  by  using  several  LCR  sections. 

For  situations  where  an  equivalent  circuit  is  needed  over  only  a  limited  range  of  frequencies. 
Reference  [9]  reported  an  alternative  to  the  use  of  a  sophisticated  software.  This  approach  is 
suitable  for  deriving  a  simple  equivalent  circuit  by  means  of  a  rational  function  representation  of 
the  input  impedance: 
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(15) 


Vis) 

Us) 


M 


i=0 

N 


where  s  =  jlTtf  (for /spanning  the  desired  range  of  frequencies),  and  the  summation  orders  Af 
and  are  selected  by  the  user  to  achieve  a  balance  between  the  accuracy  of  the  impedance-fitting 
process  and  the  complexity  of  the  resulting  equivalent  circuit  The  unknown  real  coefficients  a-, 
and  fc*  can  be  determined  by  using  either  a  least-squares  minimization  procedure  or  an 
interpolation  method  to  achieve  a  best-fit  to  the  variation  of  impedance  with  frequency  obtained 
from  the  FD-TD  modeling.  Subsequently,  these  coefficients  are  used  to  synthesize  an  equivalent 
circuit  that  is  valid  in  the  desired  range  of  frequencies. 


3.5.3  Choice  of  Methods 

Reference  [6]  provided  data  indicating  that  certain  important  distributed  components,  such  as 
ground  pads  connected  to  the  ground  plane  using  vias,  cannot  be  easily  modeled  using  the  fitting 
procedure  of  Section  3.6.2.  The  key  data  provided  (see  Fig.  8  of  [6])  indicate  that  the  phase  of 
Sii  (equivalently,  the  phase  of  Zj,)  can  have  a  dramatically  different  variation  with  frequency 
depending  upon  the  choice  of  the  reference  plane.  This  sensitivity  to  the  observation  point 
makes  it  very  hard  to  find  an  equivalent  inductance  that  is  independent  of  the  location  of  the 
observation  plane  and  also  independent  of  frequency. 

On  the  other  hand,  [6]  reported  that  the  inductance  calculated  using  the  flux/current  method 
was  found  to  be  essentially  independent  of  frequency  (see  Fig.  9  of  [6])  and  "very  accurate."  For 
such  components,  [6]  stated  that: 

As  a  result,  the  use  of  the  scattering  parameters  for  the  derivation  of  a  simple 
equivalent  circuit  becomes  much  less  successful  as  it  compares  to  the  previously 
mentioned  approach  which  is  based  on  the  flux/current  definition. 

Therefore,  the  fitting  procedure  of  Section  3.6.2  may  not  be  as  robust  as  the  more  basic 
flux/current  definition  of  inductance  of  Section  3.6.1. 

3.6  Digital  Signal  Processing  and  Spectrum  Estimation  Techniques. 

Typical  FD-TD  models  of  high-speed  circuit  structures  use  a  grid  resolution  A  that  is  dictated  by 
the  very  fine  dimensions  of  the  circuit  board  layers,  vias,  and  similar  components.  This 
resolution  is  almost  always  much  finer  than  needed  to  resolve  the  smallest  spectral  wavelength 
propagating  in  the  circuit.  As  a  result,  with  the  time  step  Ar  bound  to  A  by  numerical  stability 
considerations,  it  is  very  common  to  run  FD-TD  simulations  for  tens  of  thousands  of  time  steps 
in  order  to  fully  evolve  the  impulse  responses  needed  for  calculating  impedances,  iS-parameters, 
and  resonant  frequencies  of  passive  structures  in  electronic  circuits  operating  at  UHF  and 
microwave  frequencies. 

A  brute-force  way  to  mitigate  the  burden  of  lengthy  FD-TD  simulations  is  to  simply  truncate 
a  run  before  the  impulse  response  has  fully  evolved.  However,  this  truncation  has  the  effect  of 
viewing  the  true  time-domain  response  through  a  rectangular  window  of  duration  T  =  At. 
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In  the  frequency  domain,  this  windowing  is  translated  into  the  convolution  of  the  true  spectrum 
with  the  function  sin  (/)//.  This  convolution  widens  the  peaks  in  the  spectral  response,  causes 
other  distortions,  and  can  mask  weak  spectral  signatures.  Distortion  can  be  reduced  and  the 
spectral  resolution  increased  only  by  lengthening  the  window  duration.  Unfortunately,  this 
causes  FD-TD  computer  simulation  times  for  several  important  classes  of  problems  (particularly 
resonators)  to  be  so  lengthy  as  to  be  virtually  prohibitive. 

Recently,  there  has  been  significant  progress  in  addressing  this  problem.  The  most  fruitful 
approach  has  been  to  apply  contemporary  analysis  techniques  from  the  discipline  of  digital  signal 
processing  and  spectrum  estimation.  The  strategy  is  to  extrapolate  the  electromagnetic  field  time 
waveform  by  10:1  or  more  beyond  the  actual  TO-TD  time-stepping  window,  allowing  a  very 
good  estimate  of  the  complete  system  response  with  90%  or  greater  reduction  in  the  FD-TD 
computation  time.  In  effect,  a  low-computational-burden  extrapolation  process  at  a  limited 
number  of  reflection  or  transmission  observation  points  in  the  grid  replaces  the  computationally 
intensive  FD-TD  process  that  is  necessarily  applied  at  every  point  in  the  grid.  Then  the  reflected 
and  transmitted  field  spectra  and  the  associated  impedances  and  iS-parameters  at  the  observation 
points  can  be  efficiently  obtained  by  FFT  of  the  extrapolated  waveforms  or,  in  certain  cases, 
from  the  coefficients  of  the  extrapolation  process  itself. 

This  section  discusses  a  number  of  digital  signal  processing  and  spectrum  estimation 
techniques  that  have  appeared  in  the  recent  FD-TD  literature.  While  good  results  have  been 
obtained  with  each  of  these  methods  for  appropriate  classes  of  modeling  problems,  there  are 
emerging  tradeoff  considerations  as  FD-TD  practitioners  become  more  familiar  with  the  state  of 
the  art  of  modem  signal-processing  theoiy.  A  goal  of  this  section  is  to  illuminate  the  tradeoffs  as 
currently  understood. 

3.6.1  Pronv's  Method 

Prony's  method  [10]  has  been  used  to  extrapolate  transmission-line  matrix  (TLM)  and  FD-TD- 
computed  waveforms  for  microwave  circuits  by  a  number  of  investigators  [1 1-13].  Because  it  is 
a  technique  for  modeling  sampled  data  as  a  linear  combination  of  complex  exponentials,  it  is 
particularly  suitable  for  calculating  the  resonant  frequency  and  Q  of  a  resonating  structure,  since 
the  impulse  response  of  such  a  structure  is  characterized  by  a  superposition  of  decaying 
exponentials.  While  Prony's  method  is  not  a  spectral  estimation  technique  in  the  usual  sense,  it  is 
closely  related  to  the  least-squares  linear  prediction  methods  used  for  parametric  spectral 
estimation  [12]. 

To  begin  our  discussion,  let  us  assume  the  existence  of  N  equally  spaced  time  samples  of  the 
FD-TD-computed  impulse  response  of  a  reflected  or  transmitted  field  component  /at  observation 
point  (/,  j,  A:)  in  a  high-speed  circuit: 


J  “  Jki.k '  Juj, 


M+2 
*  ’ 


/I 


Af+A^-l 

ij.k 


(16) 


These  samples  constitute  an  observation  "window"  that  begins  at  time  step  M  and  ends  at  time 
step  M  +  N  -1.  (Note  that  {/}  is  obtained  by  decimating  the  actual  FD-TD  data  record  at 
point  (i,  /  k)  by  a  factor  of  10  or  greater.  The  FD-TD  data  are  very  oversampled  relative  to 
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what  is  needed  for  Prony's  method  and  the  other  digital  signal  processing  methods  discussed  in 
this  section.)  Now  let  each  time  sample  in  {/}  be  approximated  by  a  sum  of  p  exponentials: 


/=l 


(17) 


where  0^/i-Af<N-l,  C,  is  the  complex-valued  amplitude  of  the  fth  mode. 
At  =  {cct+ jin  ft)  At  specifies  the  modal  damping  factor  and  frequency  of  the  f  th  mode,  and  p  is 
the  order  of  the  model. 

The  direct  solution  of  (17)  is  a  difficult  nonlinear  least-squares  problem.  An  alternative 
solution  is  based  on  Prony's  method.  This  is  a  two-step  procedure  [10]  which  solves  two 
sequential  sets  of  linear  equations  with  an  intermediate  polynomial  rooting  step  that  concentrates 
the  nonlinearity  of  the  problem.  Following  [14],  the  first  step  is  to  set  up  N  —  p  equations  for 
the  p  values  of 


Jkhk 

+  A /I":.'" 

+  ...  + 

=  0 

j  Af+p+1 

+  ac/ 

+  Acr' 

=  0 

/I,-;/-  +  A 4/i'r’  = « 


W+iV-3 


iM+AT-l-p 


(18) 


This  overdetermined  system  is  solved  using  a  least-squares  algorithm  for  the  {A^}.  This  permits 
the  {//f }  to  be  found  as  the  roots  of  the  polynomial 

pf  +  Ai^^'‘  -I-  +  ...  +  Ap_iP  +  Ap  =  0  (19) 

The  second  step  is  to  find  the  {  Q}  in  (17).  This  is  accomplished  by  writing  out  (17)  for 
each  of  the  time  samples,  thereby  obtaining  a  set  of  N  equations  in  the  values  of  the  C^: 


+  Cj 

+  ...  + 

=  4". 

+ 

+  . . .  +  PpCp 

_  ^Af+1 

+  . . .  +  ^p^)  Cp 

_  -|M+2 

(20) 
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This  overdeterrained  system  is  solved  using  a  least-squares  algorithm  for  the  {Cj.  Now  all  of 
the  parameters  in  (17)  are  known,  and  it  becomes  very  easy  to: 

1.  Write  down  by  inspection  the  set  of  resonant  frequencies  and  Q-factors 
{Qt  =  \jt  ft!  directly  from  the  Prony  parameter  sets  just  obtained; 

2.  Calculate  the  complete  estimated  impulse  response  at  times  n  »  Af  +  N  - 1 
until  it  is  observed  to  decay  to  zero,  at  which  point  an  FFT  provides  the 
magnitude  and  phase  of  the  frequency-domain  transfer  function  with  no 
windowing  artifacts. 

Reference  [13]  reported  results  using  Prony’s  method  for  the  resonant  frequencies  and  Q’s  of 
the  three  lowest  TEo„  modes  of  a  cylindrical  dielectric  resonator  of  permittivity  e,.  =  38,  radius 
a  =  5.25  mm,  and  length  I  =  4.6  mm.  A  cylindrical  FD-TD  grid  was  used,  with  parameters 
Ar  =  0.29167  mm,  Az  =  0.2875  mm,  and  At  =  0.4795  ps.  A  time-stepping  window  extending 
over  about  2*’  (131,072)  steps  was  needed  to  obtain  converged  results  for  the  Q  of  the 
mode  when  using  an  FFT  applied  to  the  windowed  FD-TD  data.  However,  using  Prony’s  method 
of  order  p  =  10  and  a  decimation  factor  of  50,  only  3000  time  steps  were  required  to  obtain  the 
same  Q,  a  reduction  of  about  98%.  In  fact,  reductions  can  be  even  greater  in  the  case  of 
structures  having  even  higher  Q  or  with  very  closely  spaced  resonant  frequencies. 

Prony’s  method  shows  good  resolution  with  relatively  short  data  sequences  and  presents  no 
windowing  problem.  The  main  difficulty  with  this  approach  lies  in  its  determination  of  p,  the 
user-selected  order  of  the  model.  If  the  value  of  p  is  less  than  the  number  of  actual  modes 
excited  in  the  structure,  the  spectral  resolution  is  poor.  However,  if  p  is  selected  to  be  too  high, 
spurious  (nonphysical)  modes  appear.  The  following  are  two  guidelines  given  in  [13]  to  deal 
with  these  problems: 

•  Before  Prony's  method  is  applied,  digitally  low-pass  filter  the  time-domain  response  to 
limit  the  number  of  resonant  modes  and  therefore  the  number  of  parameters  to  calculate. 

•  If  spurious  modes  are  suspected,  apply  Prony’s  method  with  the  time  sequence  of  the 
sample  points  in  reverse  order.  Here  the  real  modes  appear  with  positive  damping  factors 
(fitf  >  0),  but  the  spurious  modes  have  negative  damping  factors  (a^  <  0)  [15]. 

3.6.2  AwtoregressivsJdDdeis 

The  most  popular  time  series  modeling  approach  used  in  modem  spectral  estimation  is  the  class 
of  linear  predictors  or  autoregressive  (AR)  models.  This  is  because  an  accurate  estimate  of  the 
AR  parameters  can  be  derived  by  solving  a  set  of  linear  equations.  In  contrast  to  Prony’s  method, 
which  uses  a  sum  of  deterministic  exponential  functions  to  fit  the  data,  the  AR  approach 
constructs  a  random  model  to  Ht  a  statistical  data  base  to  the  second-order. 

We  again  assume  the  FD-TD-computed  impulse  response  of  (16),  i.e.,  N  equally  spaced  time 
samples  beginning  at  time  step  M.  This  time  series  is  said  to  represent  the  realization  of  an  AR 
process  of  order  pHii  satisfies  the  following  relation: 


=  -o./!”;;  -  ‘hG 


ii-2 


(21) 
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where  the  constants  Oj ,  Oj , . . are  the  AR  parameters  to  be  determined  from  {/},  and  q{n) 
is  a  white  noise  process  whose  variance  also  has  to  be  found  to  carry  out  the  extrapolation  of 
{/}.  It  is  clear  that  (21)  permits  the  prediction  of  a  new  value  of  the  time  series  from  known 
previous  values  once  the  AR  parameters  are  in  hand. 

There  exists  many  published  approaches  to  evaluating  the  a,  [16].  In  fact,  a  primary 
categorization  of  AR  techniques  is  by  the  means  used  to  calculate  these  parameters.  In  this 
section,  we  consider  three  specific  recent  applications  of  AR  to  the  extrapolation  of  FD-TD 
electromagnetic  field  records,  categorized  according  to  this  criterion:  (1)  the  covariance  method, 
(2)  the  forward-backward  method,  and  (3)  the  nonlinear  predictor.  The  choice  of  approach  is 
very  important,  since  these  three  methods  have  substantially  different  operational  characteristics. 

Covariance  Method 

Reference  [17]  reported  the  use  of  the  covariance  method  for  determining  the  AR  parameters. 
This  method  involves  setting  up  and  solving  the  following  py.  p  linear  system  of  equations: 


•c^/1,1)  c,,(l,2) 

c,/2,l)  c,/2,2) 

Cff{pA)  Cjf{p,2) 

where  are  covariances  defined  by 

C  ff  (fiCf  P)  — 


Cffilp)' 

■c^/(l,0)- 

Cff(2,  p) 

«2 

=  - 

0) 

1 _ 

1 

o 

1 _ 

(22) 


(23) 


The  matrix  of  (22)  is  Hermitian  and  positive  semidefinite,  and  can  be  solved  by  Cholesky 
decomposition. 

As  stated  in  [17],  a  key  issue  that  arises  when  using  the  covariance-based  AR  model  is 
choosing  the  order  of  the  model,  p.  The  use  of  a  low-order  AR  model  of  this  type  causes  the 
extrapolated  waveform  to  attenuate  quickly  in  a  nonphysical  manner.  To  compensate,  a  high- 
order  model  is  required.  However,  a  high-order  model  can  cause  divergence  problems  in  some 
cases  because  of  statistical  instabilities  introduced  by  the  large  order.  These  competing  effects 
cause  the  results  of  the  covariance-based  AR  model  to  be  sensitive  to  the  model  order, 
diminishing  its  usefulness  somewhat.  Two  commonly  used  ways  to  estimate  p  are  the  final 
prediction  error  technique  and  the  Akaike  Information  Criterion  [18].  Both  are  based  upon  the 
estimated  predictor  error  power  and  are  regarded  as  general  guides  for  AR  model  selection. 

Using  the  covariance-based  AR  method  of  (22)  and  (23),  [17]  reported  the  successful 
extrapolation  to  over  30,000  time  steps  of  FD-TD  field  records  spanning  only  2000  to  3500  time 
steps.  A  decimation  factor  of  10  was  used,  and  the  order  of  the  AR  model  was  50.  Another 
example  discussed  in  [17]  involved  modeling  a  strip  line  double-stub  structure  that  required  a 
substantial  computation  space  because  of  the  large  ratio  of  the  maximum  to  minimum 
dimensions  of  the  structure.  Here,  the  lengthy  wave  propagation  delays  along  the  stubs  required 
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extrapolating  the  FD-TD  data  to  100,000  time  steps  to  achieve  convergence  of  the  transmission 
coefficient.  This  extrapolation  was  accomplished  using  a  91st-order  AR  model  working  from  an 
FD-TD  data  base  covering  iterations  3 100  to  8900,  and  decimated  by  a  factor  of  25.  Although 
very  good  agreement  was  observed  with  benchmark  frequency-domain  data,  the  results  of  the 
AR  model  were  found  to  be  very  sensitive  to  the  order  of  the  model,  much  more  so  than  the 
50th-order  model  used  for  the  previous  two  examples. 


Forward-Backward  Method 


The  results  of  a  covariance-based  AR  model  can  be  sensitive  to  the  order  of  the  model,  because 
its  parameters  often  have  low  accuracy.  The  low  accuracy  results,  in  turn,  because  the 
ensemble  calculation  that  is  rigorously  required  to  find  the  covariance  is  instead  substituted  in 
these  methods  by  using  the  law  of  large  numbers  and  by  approximating  the  covariance  with 
inexact  functions  of  the  known  time  signal.  A  useful  alternative  approach  based  on  an 
unconstrained  least-squares  estimation  of  the  AR  parameters  was  proposed  independently  in  [19] 
and  [20]  and  calculated  more  efficiently  in  [21].  This  approach,  called  forward  and 
backward  prediction  method  [20],  was  used  for  the  first  time  in  [22]  to  extrapolate  FD-TD- 
computed  field  records.  It  avoids  the  problems  of  the  covariance-based  AR  models  by  working 
directly  with  the  time-domain  data,  rather  than  calculating  the  covariance  functions  of  the  data. 

Reference  [21]  reported  setting  up  the  following  (p  + 1)  x  (p  + 1)  linear  system  of  equations 
for  the  forward-backward  method: 


>(0,0)  ...  r(0,p)' 

‘l  * 

«i 

0 

•  •  • 

•  •  • 

— 

0 

0 

r(p,  0)  ...  r(p,  p) 

0 

where,  for  0  ^  Cf,  ^  p,  we  have 


(24) 


(25a) 


(25b) 


Reference  [21]  proved  that  matrix  [r(a,  p)\  of  (24)  is  always  positive  semidefinite;  and  if  it  is 
nonsingular,  then  it  is  positive  definite.  This  is  the  most  common  case,  since  otherwise  the 
recursive  determination  for  the  AR  parameters  would  not  be  possible.  Reference  [21]  also 
derived  an  elegant  algorithm  to  solve  the  linear  system  of  (24)  without  inverting  [r(a,  j5)], 
thereby  requiring  only  0{p^)  arithmetic  operations.  However,  this  sophistication  does  not 
appear  to  be  necessary  when  extrapolating  just  a  few  FD-TD  records  having  a  size  (after 
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decimation)  of  perhaps  a  few  hundred  data  points,  and  an  AR  model  order  p  of  50  or  less.  In 
fact,  [19, 20]  successfully  solved  (24)  for  the  {a}  using  simple  matrix  inversion. 

Reference  [21]  reported  favorable  results  for  the  forward-backward  method  versus  existing 
popular  AR  approaches  such  as  the  Burg  algorithm  and  the  Yule-Walker  algorithm.  The 
generated  AR  spectra  had  less  bias  in  the  frequency  estimate  of  the  spectral  components,  reduced 
variance  in  frequency  estimates  over  an  ensemble  of  spectra,  and  no  spectral  line  splitting.  Much 
of  the  performance  improvement  was  attributed  in  [21]  by  the  removal  in  the  forward-backward 
method  of  the  constraint  that  the  estimation  procedure  always  place  poles  inside  the  unit  circle  in 
the  z-transform  plane.  It  is  clear  that  the  removal  of  this  constraint  presents  a  potential 
instability,  indicating  the  need  for  caution  when  using  this  method.  However,  based  upon  the 
results  of  over  3000  spectra  with  varying  signal-to-noise  ratios  generated  from  actual  data 
records  via  the  forward-backward  method,  [21]  noted  that  fewer  than  1%  had  poles  outside  of 
the  unit  circle.  In  addition,  [20]  found  no  unstable  poles  when  simulating  narrowband  processes. 
Overall,  these  findings  indicate  that  the  forward-backward  method  is  sufficiently  robust  for 
practical  engineering  application  and  is  an  attractive  alternative  to  the  Burg  and  related 
algorithms  for  AR  spectral  estimation. 

Reference  [22]  was  the  first  to  report  the  application  of  the  forward-backward  method  to  AR 
extrapolation  of  FD-TD  data  records  and  corresponding  spectral  estimation.  Results  were  given 
for  a  test  problem  involving  a  rectangular  PEC  plate  embedded  in  a  high-permittivity  half  space 
and  excited  by  a  sinusoidal  point  source  located  in  the  free-space  region.  Both  an 
autocorrelation-based  Yule-Walker  AR  model  and  the  forward-backward  AR  model  were 
"trained"  on  the  same  FD-TD  field  waveform  record,  data  points  from  201  to  800  time  steps 
("appropriately  decimated  and  low-pass  filtered").  Both  models  were  then  used  to  generate  an 
extrapolated  sequence  extending  over  5000  time  steps.  The  forward-backward  method  exhibited 
two  key  advantages  relative  to  the  Yule-WaUcer  method:  (1)  It  provided  more  accurate  spectra, 
and  (2)  Its  order  was  much  lower,  p  =  5  versus  p  =  50.  It  was  concluded  in  [22]  that  the 
forward-backward  AR  model  appears  to  provide  a  very  good  approach  to  FD-TD  signal 
extrapolation  and  spectrum  estimation  both  in  terms  of  accuracy  and  speed. 

Nonlinear  Predictor 

In  addition  to  the  use  of  covariance-based  AR  linear  prediction  to  extrapolate  FD-TD  waveform 
records  (discussed  earlier),  [17]  reported  the  application  of  a  nonlinear  predictor  for  the  same 
purpose.  The  basic  form  of  this  predictor  was  the  exponential  autoregressive  (EXPAR)  model  of 
order  p  given  by  [23,24]: 

f\lhk  =  +  error  term  (26) 

where  {a},  {fc),  and  S  are  the  modeling  parameters  that  are  "trained"  from  the  windowed  and 
decimated  FD-TD  field- versus-time  record. 

The  EXPAR  model  was  applied  in  [17]  to  the  same  large  double-stub  strip  line  problem 
discussed  earlier  in  [17],  and  similarly  used  to  extrapolate  an  8900-step  FD-TD  data  record  to 
100,000  time  steps.  The  nonlinear  EXPAR  model  provided  excellent  accuracy  for  S21  over  a 


wide  dynamic  range  upon  Fourier  transformation  of  the  100,000-step  record.  Moreover,  the 
EXPAR  results  were  much  less  sensitive  to  the  model  order  than  were  the  results  of  the  91st- 
order  covariance-based  linear  AR  model.  However,  the  reported  order  of  the  EXPAR  model  in 
this  paper  exceeded  80,  still  very  high.  A  necessary  line  of  inquiry  involves  applying  the 
forward-backward  linear  AR  model  to  problems  of  this  size  and  complexity  to  establish  whether 
its  superior  characteristics  can  markedly  reduce  the  required  model  order. 

3.6.3  System  Identification 

For  a  linear  system,  it  is  well  known  that  the  output  signal  y(f)  is  related  to  the  input  signal  x(0 
by  the  convolution 


y(t)  =  h(t)  *  xit)  (27a) 

where  h(t)  is  the  system  impulse  response.  Rigorously,  convolution  requires  knowledge  of  the 
entire  time  history  of  xit).  For  a  discrete-time  linear  system,  the  time-domain  output/input 
relation  of  (27a)  can  be  expressed  as  the  following  nonrecursive  weighted  sum  of  time  samples 
of  xit),  assuming  zero  initial  conditions: 


,>1”  = 

This  operation  typically  characterizes  what  is  called  a  finite-duration  impulse  response  (FIR) 
fdter. 

Now  consider  an  alternate  system  identification  (SI)  approach  [25-28]  to  characterize  a 
discrete-time  linear  system.  This  approach  specifies  the  time-domain  output/input  relation  as  the 
following  recursive  difference  equation: 


(28) 

where  the  system  parameter  sets  {a}  and  {b}  are  initially  unknown  and  must  be  "trained"  on  a 
windowed  and  decimated  portion  of  the  input  signal  record.  Because  the  SI  approach  has 
feedback,  it  is  termed  an  infinite-duration  impulse  response  (HR)  filter.  Reference  [28]  notes 
that  it  is  known  that  the  HR  filter  can  provide  a  less  complex,  lower-order  realization  than  the 
FIR  filter  for  low-pass  systems.  Therefore,  the  SI  HR  formulation  of  (28)  should  be  more 
efficient  for  processing  FD-TD  data  than  the  convolutional  FIR  formulation  of  (27b),  since 
FD-TD  is  inherently  a  low-pass  system  (i.e.,  its  spectrum  of  space/time  Fourier  modes  is 
bounded  by  the  Nyquist  limit  established  by  the  grid  sampling). 

We  now  consider  the  search  procedure  in  ({a},  {ft})  space  reported  by  [27]  to  obtain  the 
optimum  system  parameter  sets.  In  the  context  of  FD-TD,  the  process  deals  with 
electromagnetic  field  quantities,  so  we  use  the  field  component  f  instead  of  y,  and  the  usual 
space-time  notation.  Equation  (28)  is  now  written  in  compact  form  as  the  dot  product  of  two 
one-dimensional  arrays: 
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(29) 


/Hm  =  Kllf  •  hL,.] 

Here  [<I>]  is  a  data  array  containing  the  past  and  present  values  of  the  input  and  output  field 
samples  at  grid  point  (i,  y,  k),  [Gg]  contains  the  parameter  sets  {a}  and  {b}  unique  to  this  grid 
point,  and  T  denotes  the  array  transpose  operation.  However,  at  the  beginning  of  the  search 
procedure  (say,  time  step  n),  the  elements  of  [Gq]  are  unknown.  We  make  the  arbitrary  guess 
[Gfl].  (Reference  [27]  reported  success  with  the  initial  guess  [Qq]  =  0,  that  is,  the  search 
beginning  at  the  origin  of  the  ({a},  {fc})  space.)  ^Using  the  existing  input  and  output  data  records 
at  time  step  n,  an  estimate  for  the  output  field  n  can  now  be  obtained  in  terms  of  the  guessed 
system  parameters: 


K»r-[C] 


(30) 


The  difference  between  the  actual  and  estimated  FD-TD  field  outputs  at  time  step  n  is  obtained 
by  subtracting  (30)  from  (29),  and  then  minimized  with  respect  to  the  elements  of  {a}  and  {b}. 
This  yields  the  following  system  parameter  update  relations: 
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(32) 


where  P|.^  ^  =  /,  the  identity  matrix.  The  procedure  of  (31)  and  (32)  is  continued  until 
discrepancy  between  the  estimated  output  and  the  FD-TD-computed  field  value,  is  sufficiently 
small.  At  this  point,  the  system  parameters  at  point  (j,  j,  k)  have  been  "trained,"  and  a  useful 
extrapolation  of  the  field  record  is  possible.  From  a  mathematics  perspective,  this  method 
provides  an  orthogonal  projection  search  in  the  system  parameter  space  which  usually  results  in 
rapid  parameter  convergence.  It  is  properly  called  a  "least-squares-based  system  identification 
projection  algorithm  for  a  deterministic  autoregressive  moving  average  (ARMA)  model." 

An  ARMA  model  with  system  parameters  N  =  40  and  N'  =  40  in  (28)  was  applied  in  [27] 
to  calculate  the  resonant  frequencies  of  a  PEC  rectangular  cavity  partially  filled  with  a  dielectric 
slab.  The  cavity  was  excited  at  its  center  plane  by  imposing  a  pulsed  TE,o  mode  field 
distribution,  and  the  time  waveform  of  this  excitation  was  used  as  the  input  signal  x(t)  in  (28). 
The  parameter  search  began  at  n  =  0  with  {a}  =  {0},  and  the  parameter  values  were  updated  at 
each  FD-TD  time  step  (no  decimation).  Upon  convergence  of  these  parameters,  the  resonant 
frequencies  of  the  cavity  could  be  derived  with  no  Fourier  transform  of  the  output  signal,  if 
desired.  Instead,  the  2^transform  of  (28)  would  be  evaluated  on  the  unit  circle  which  is  defined 
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in  terras  of  the  system  parameters,  thereby  obtaining  the  resonant  frequencies  directly  from  the 
poles  of  the  ARMA  model. 

Overall,  the  ARMA  SI  algorithm  requires  minimal  additional  computational  cost.  It 
provides  a  rapidly  converging  recursive  processing  of  the  input  and  output  values  of  the  FD-TD 
simulation,  and  can  be  performed  concurrently  with  the  usual  FD-TD  time-stepping  at  any 
number  of  observation  points  in  the  grid.  In  addition  to  permitting  an  extrapolation  of  the  local 
field  values  at  the  observation  points,  the  frequency  spectrum  is  immediately  available  without 
needing  to  apply  a  Fourier  transform.  Further,  the  input/output  nature  of  the  SI  formulation 
permits  an  efficient  natural  treatment  of  diakoptic  spUtting  of  the  computation  space  without 
needing  to  use  convolutions  [28, 29]. 

4.0  TASK  2.  PFRFORM  VTSTTALTZATIONS  RESULTING  FROM  FD-TD  MQDBLMS 

This  section  reports  selected  visualizations  of  propagating  signals  within  commercial  circuit 
boards  of  contemporary  design.’  Because  we  are  working  in  the  time  domain,  we  can  obse^e 
the  dynamics  of  potentially  important  physical  phenomena  such  as  reflections,  non-uniformities 
of  current  flow,  edge  and  comer  enhancements,  and  fields  propagating  through  the  dielectric 
media  rather  than  flowing  along  the  desired  metallic  signal  path. 

4.1  Example  1:  Coupling  Between  ParaUel  Adjacent  Via  HoS 

This  example  involves  the  detailed  3-D  FD-TD  modeling  of  a  very  high-speed  computer  circuit 
board  having  more  than  twenty  0.004-inch  thick  layers.  The  board  is  penetrated  by  a  regular 
array  of  0.012-inch  diameter  circular  via  pins  located  on  0.1-inch  centers.  A  uniform  grid  cell 
size  of  A  =  0.004  inch  is  used  in  the  FD-TD  model,  with  a  time  step  Ar  =  0.169  ps.  'Die  power 
and  signal  metallizations  of  the  circuit  board  are  assumed  to  be  monolayers  of  tangential  electric 
field  components  set  to  zero. 

Fig.  1  is  a  color  visualization  of  the  early-time  response  of  the  circuit  board  to  a  0.1 -ns 
Gaussian  current  pulse  propagating  down  a  single  vertied  via  pin.  The  via  pin  (in  a  circular  via 
hole)  is  excited  by  pulsing  a  vertical  electric  field  component  in  the  FD-TD  grid  just  above 
the  pin  and  below  a  simulated  ground  strap  connected  to  the  ground  return.  The  excitog  pulse  is 
selected  to  have  a  bandwidth  comparable  to  that  of  the  actual  stream  of  nanosecond  digital  pulses 
running  through  the  via  (including  the  harmonic  energy  in  these  pulses).  Fig.  1  depicts  the 
strength  of  the  magnetic  field  in  a  vertical  cut  through  the  board,  using  a  color  scale  of  red  - 
yellow  -  green  -  blue  -violet  to  visualize  maximum  to  minimum  field  intensities,  respectively. 
We  observe  substantial  pin-to-pin  crosstalk  due  to  strong  coupling  of  the  magnetic  field  from  the 
excited  via  pin  to  the  adjacent  unexcited  via  pins.  In  fact,  the  peak  coupled  magnetic  field  (and 
thus  the  peak  coupled  pin  current)  is  down  by  only  about  2  to  3  dB  in  Ae  via  pins  immediately 
adjacent  to  the  excited  pin.  This  represents  a  potentially  very  troublesome  "self  jamming" 
problem  inherent  in  this  circuit  board  design. 


*The  names  of  the  companies  involved  arc  not  relevant  to  this  discussion  and  are  therefore  omitted. 
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Fig.  1  Color  visualizalion  of  the  early-time  response  of  a  >20-Iaycr  circuit  board  to  a  0.1-ns  Gaussian  current 
pulse  propagating  down  a  single  vertical  via  pin.  The  figure  depicts  the  strength  of  the  magnetic  field  in  a 
vertical  cut  tlirough  the  board,  using  a  color  scale  of  red  -  yellow  -  green  -  blue  -violet  to  visualize 
maximum  to  minimum  field  intensities.  ITiere  is  substantial  pin-to-pin  crosstalk  due  to  strong  coupling  of 
the  magnetic  field  from  the  excited  center  via  pin  to  the  two  adjacent  unexcited  via  pins. 
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4.2  F.xample  1:  Coupling  Between  Perpendicular  Circuit  Traces 

This  example,  constructed  especially  for  this  Rome  Laboratory  program,  involves  the  detailed 
3-D  FD-TD  modeling  of  coupling  from  a  single  signal  trace  in  a  commercial  very  high-speed 
computer  circuit  board  to  three  perpendicular  signal  traces  in  the  adjacent  metalization  layer.  A 
uniform  grid  cell  size  of  A  =  12.7  |xm  (equal  to  0.0(X)5  inch,  or  0.5  mil)  is  used  in  the  FD-TD 
model  with  a  time  step  Ar  =  20  fs.  As  in  Example  1,  all  signal  traces  of  the  circuit  board  are 
assumed  to  be  monolayers  of  tangential  electric  field  components  set  to  zero. 

Figs.  2a  and  2b  depict,  respectively,  the  cross-section  and  plan  views  of  this  geometry.  The 
yellow  region  in  each  figure  denotes  the  interior  of  the  modeled  circuit  board  filled  with  a 
dielectric  of  er=  3.0,  and  the  blue  lines  denote  metallizations.  The  array  of  dots  represents  grid 
reference  points  spaced  at  5-cell  (63.5  ^im  or  2.5  mils)  increments.  There  is  nothing  physically 
located  at  these  dots. 

Referring  to  Fig.  2a,  an  approximately  30-ps  duration  Gaussian  signal  pulse  is  assumed  to  be 
incident  into  the  page  (in  the  +z  direction)  on  the  lower  trace,  which  is  3  cells  (38.1  pm  or  1.5 
mils)  wide  and  2  cells  (25.4  p.m  or  1  mil)  above  the  lower  ground  plane.  There  are  three 
crossing  traces  located  in  a  plane  2  cells  (25.4  pm  or  1  mil)  above  the  signal  trace  and  oriented 
perpendicular  to  it  in  the  x-direction.  The  upper  ground  plane  is  located  2  cells  (25.4  pm  or 
1  mil)  above  the  plane  of  the  three  crossing  traces. 

Referring  to  Fig.  2b,  each  of  the  crossing  traces  is  also  3  cells  (38.1  pm  or  1.5  mils)  wide, 
with  a  center-to-center  spacing  of  25  cells  (317.5  pm  or  12.5  mils)  in  the  z-direction  relative  to 
the  adjacent  crossing  trace.  The  first  of  these  crosses  the  signal  trace  25  cells  (317.5  pm  or  12.5 
rails)  from  the  source.  All  traces  are  terminated  in  a  4-ceU  thick  Berenger  perfectly  matched 
layer  (PML)  absorbing  boundary  condition  (ABC). 

Figs.  3a  and  3b  visualize  a  time  sequence  of  snapshots  of  the  longitudinal  current  density 
distributions  on  the  signal  and  crossing  traces  as  indicated  by  the  transverse  tangential  H-field 
component  at  the  surface  of  each  trace.  (The  observation  perspective  is  the  plan  view  of  Fig.  2b.) 
The  upper  panels  show  Hx  at  the  plane  of  the  signal  trace,  thereby  providing  Jz  along  the  trace. 
The  lower  panels  show  the  corresponding  Hz  at  the  plane  of  the  crossing  traces,  thereby 
providing  Jx  along  these  traces.  The  color  scale  provides  bipolar  information,  with  blue  - 
light  blue  -  white  denoting  increasingly  positive  values,  and  red  -  yellow  -  white  denoting 
increasingly  negative  values.  We  see  that  the  passage  of  the  +z-directed  incident  pulse  on  the 
signal  trace  successively  induces  a  rightward  (+x-directed)  pulse  and  a  leftward  (-x-directed) 
pulse  on  each  of  the  crossing  traces  above  it  which  propagate  away  from  the  vicinity  of  the  signal 
trace.  The  resulting  current  density  values  on  each  crossing  trace  have  odd  symmetry  in  the 
x-direction  about  a  point  just  above  the  signal  trace.  As  the  signal  pulse  leaves  the  grid,  so  do 
the  coupled  pulses  on  the  crossing  lines,  due  to  the  action  of  the  PML  ABC. 

Fig.  4  graphs  the  spatial  distribution  of  Hx  (equivalently,  along  the  signal  trace  (Line  1) 
and  Hz  (equivalently,  /,)  along  the  first  crossing  trace  (Line  2)  at  approximately  the  time  of 
maximum  coupling.  Fig.  5  graphs  the  corresponding  spatial  distributions  of  the  voltage  of  each 
line,  referenced  to  the  lower  ground  plane.  The  peaJc  coupling  from  the  signal  trace  to  the  first 
crossing  trace  is  about  -28  dB.  This  is  considerably  weaker  than  the  -2  or  -3  dB  coupling  noted 
in  Example  1  (adjacent  parallel  via  pins),  and  should  not  present  an  electromagnetic 
compatib^ty  problem  even  at  contemporary  supercomputer  clocks  of  about  1  GHz. 
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Fig.  2a  Vertical  cross-section  of  3-D  FD-TD  model  of  coupling  from  a  single  z-directed  signal  trace  in  a 
commercial  very  high-speed  computer  circuit  board  to  three  jc-directed  signal  traces.  Yellow  denotes  the 
interior  of  the  circuit  board  filled  with  a  dielectric  of  = 3 .0;  the  blue  lines  denote  metallizations.  The  dots 
are  grid  reference  points  spaced  at  5-cell  (63.5  |im  or  2.5  mils)  increments;  there  is  nothing  physically 
located  at  these  dots.  Vertical  separation  between  each  metalization  plane  =  2  cells  (25.4  |jm  or  1  mil). 
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Fig.  2b  Plan  view  of  3-D  FD-TD  model  of  coupling  from  a  single  z-directed  signal  trace  in  a  commercial  very 
high-speed  computer  circuit  board  to  three  jr-directed  signal  traces.  Yellow  denotes  the  interior  of  the 
circuit  board  filled  with  a  dielectric  of  er=  3.0;  the  blue  lines  denote  metal  traces.  The  dots  are  grid 
reference  points  spaced  at  5-cell  (63.5  pm  or  2.5  mUs)  increments;  there  is  nothing  physically  located  at 
these  dots.  Lateral  (z)  separation  between  each  crossing  trace  =  25  cells  (317.5  pm  or  12.5  mils). 
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Fig.  3a  Time  sequence  visualizing  snapshots  of  the  longitudinal  current  densities  on  the  signal  and  crossing  traces. 
Upper  panels  show  at  the  plane  of  the  signal  trace,  yielding  along  the  trace.  Lower  panels  show  the 
corresponding  at  the  plane  of  the  crossing  traces,  yielding  along  these  traces.  Color  scale:  blue  - 
light  blue  -  white  show  increasingly  positive  values;  red  -  yellow  -  white  show  increasingly  negative  values. 
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Fig.  3b  Continuation  of  the  time  sequence  of  Fig.  3a. 
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Fig.  4  Cuirent  density  distributions  at  the  time  of  maximum  coupling  from  the  signal  trace  (Line  1)  to  the  first 
crossing  trace  (Line  2).  Top:  variation  of  //,  (equivalently,  7j)  with  z  along  the  signal  trace.  Bottom, 
variation  of  H*  (equivalenfly,  with  x  along  the  first  crossing  trace,  where  x  =  0  is  just  above  the  center 
line  of  the  signal  trace. 
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Potential  distributions  at  the  time  of  maximum  coupling  from  the  signal  trace  (Line  1)  to  the  first  trossmg 
trace  (Line  2).  Top:  variation  of  potential  with  z  along  the  signal  trace.  Bottom:  variation  of  potential  wim 
X  along  the  first  crossing  trace,  where  x  =  0  is  just  above  the  center  line  of  tbe  signal  trace.  All  potentials 
are  referenced  to  the  lower  ground  plane. 
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This  example  involves  the  detailed  3-D  FD-TD  modeling  of  an  837-MHz  resonant  RF  power 
combiner  constructed  of  a  single-layer  ceramic  circuit  board.  Field  reports  indicated  that  this 
device  is  susceptible  to  arcing  and  scorching,  but  this  failure  mode  is  sporadic  and  unpredictable. 

The  power  combiner  consists  of  an  S-shaped  stripline  laid  down  on  a  1.524-mm  (60-mils) 
thick  alumina  substrate  over  ground  having  the  dielectric  parameters  £r  ~  9-8  and  tanS  —  0.0001. 
The  combiner  has  six  straight  stripline  sections  laid  out  at  right  angles.  Stripline  lengths  are  as 
follows: 

Trace  1.  36.068  mm  (1420  mils)  z-directed;  connected  to  a  50-ohm  RF  power  source. 

Trace  2.  8.128  mm  (320  mils)  x-directed 

Traces.  19.304mm  (760 mils)  z-directed;  designated  the  "center  trace." 

Trace  4.  8.128  mm  (320  mils)  x-directed;  designated  the  "top  trace." 

Trace  5.  20.828  ram  (820  mils)  z  directed 

Trace  6.  13.208  mm  (520  mils)  x-directed;  ends  with  a  50-ohm  load  to  ground. 

To  adjust  the  operation  of  the  combiner,  a  15.24  x  20.32  x  2.032-mm  (600  x  800  x  80-mils) 

barium  titanate  "puck"  having  the  dielectric  parameters  £r  —  37  and  tan5  =  0.0001  rests  on  the 
surface  of  the  circuit  board  so  that  it  partially  covers  the  S-shape.  A  field  technician  adjusts  the 
position  of  this  ceramic  puck  to  vary  the  partial  overlap,  thereby  tuning  the  resonance  of  the 
combiner.  It  was  found  that  arcing  is  most  likely  when  the  left  edge  of  the  puck  is  located  0.508 
mm  (20  mils)  to  the  right  of  the  center  trace,  and  the  top  edge  of  the  puck  is  located  2.032  mm 
(80  rails)  beyond  the  top  trace.  By  observing  the  scorch  marks  on  the  surface  of  the  alumina 
substrate,  it  was  determined  that  the  arcing  originates  at  the  junction  of  the  center  and  top  traces 
(Traces  3  and  4,  respectively)  and  burrows  under  the  barium  titanate  puck  in  the  +x-direction 
towards  Trace  1. 

The  FD-TD  model  uses  a  uniform  space  increment  of  0.508  ram  (20  mils)  and  a  time  step  of 
0.8467  ps.  Fig.  6  is  a  color  visualization  of  the  FD-TD-computed  amplitude  of  the  sinusoidal 
steady-state  Ex  field  distribution  at  the  surface  of  the  alumina  substrate.  Ex  was  selected  for 
visualization  since  the  arcing  was  observed  to  occur  in  the  x-direction.  Note  that  the  position  of 
the  barium  titanate  puck  is  visible  as  a  weak  rectangular  shadow  in  this  figure.  From  the  color 
coding,  where  white  denotes  maximum  field  amplitude,  we  see  that  the  FD-TD  model  predicts 
that  the  maximum  Ex  field  arises  at  the  junction  of  the  center  and  top  traces,  exactly  where  arcing 
was  experienced  in  practice.  In  fact,  elevated  Ex  fields  are  seen  to  arise  all  along  the  gap  between 
the  center  trace  and  the  left  side  of  the  puck.  Further,  scaling  of  the  computed  peak  Ex 
intensities  to  the  normal  operating  power  of  the  combiner  shows  electric  field  levels  nearly  at  the 
breakdown  threshold  of  air  at  standard  temperature,  pressure,  and  humidity.  This  helps  to 
explain  the  erratic  nature  of  the  observed  arcing,  since: 

1.  The  combiners  were  not  hermetically  sealed,  and  the  variability  of  normal  atmospheric 
conditions  could  shift  the  breakdown  threshold  of  the  air  inside  the  combiner. 

2.  Slight  differences  in  the  manual  adjustment  of  the  puck  position  by  various  field 
technicians  at  various  sites  could  slightly  increase  or  decrease  the  peak  electric  field, 
thereby  causing  the  combiner  to  be  respectively  closer  to,  or  further  from,  breakdown. 


Fig.  6  Visualization  of  the  FD-TD-computed  amplitude  of  the  sinusoidal  steady-state  field  distribution  at  the 
surface  of  the  alumina  substrate  of  the  RF  power  combiner  at  837  MHz.  Color  scale:  decreasing  intensity 
from  white  to  yellow  to  red  to  purple  to  black.  Note  that  the  position  of  the  barium  titanate  puck  is  visible 
as  a  weak  rectangular  shadow  having  its  left  edge  just  to  the  right  of  the  center  trace. 
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5.0  TASK  3.  PROCEED  WITH  FD-TD/SPICE  SOFTWARE  DEVELOPMENT 


5.1  Introduction 

This  section  reports  the  continuation  of  the  development  of  a  hybrid  software  combining  FD-TD 
field  analysis  and  SPICE  nonlinear  circuit  analysis.  This  work  should  eventually  permit  full- 
wave  simulation  of  high-speed  digital  modules,  including  both  the  operation  of  the  active  logic 
devices  and  the  passive  interconnects  between  these  devices. 

Reference  [30]  demonstrated  that  it  is  possible  to  construct  local  software  links  between  the 
three-dimensional  FD-TD  Maxwell's  equations  software  and  appropriate  SPICE  kernels.  This 
software  forms  an  interface  between  the  time-stepping  of  the  linear,  distributed,  full-vector 
electromagnetic  {E,  H)  fields  and  the  nonlinear,  lumped,  scalar  voltages  and  currents  (V,  /)  of 
arbitrarily  complex  electronic  circuits  (such  as  gate  arrays)  that  have  negligible  electromagnetic 
wave  propagation  effects  within  the  circuit.  Effectively,  the  SPICE  link  permits  FD-TD  subgrid 
models  of  transistors,  logic  gates,  and  gate  arrays  to  be  constructed  that  incorporate  all  important 
aspects  of  their  circuit  physics,  including  nonlinearides  at  inputs  and  outputs,  device  parasitics, 
and  time  delays  and/or  input/output  logic-state  shifts  due  to  analog  or  digital  operations 
occurring  within  the  device. 

5.2  Formulation 

The  connection  between  Maxwell's  equations  and  SPICE  is  simple  and  elegant.  The  approach  is 
suitable  for  modeling  an  arbitrary  two-port  active  network  in  the  context  of  the  time-varying 
electromagnetic  field.  With  the  hybrid  FD-TD/SPICE  method,  the  network  is  replaced  with 
equivalent  current  sources  at  the  location  of  its  two  ports  in  the  space  grid.  At  each  port,  the 
equivalent  sources  must  satisfy  the  possibly  nonlinear  current- voltage  relationship  of  the  network 
and  provide  the  same  response  as  the  network  when  an  electromagnetic  wave  arrives. 

We  now  summarize  the  basis  of  defining  the  equivalent  sources.  Consider  the  differential 
form  of  Ampere's  Law  written  at  the  location  of  a  single,  pointwise-connected  circuit  element: 

e^+J{E)  =  VxH  (33) 

at 

Here,  the  first  term  on  the  left-hand  side  is  the  local  displacement  current  density.  The  second 
term  on  the  left-hand  side  is  the  conduction  current  density  due  to  the  circuit  element.  In  general, 
the  conduction  current  density  can  have  a  nonlinear  dependence  upon  the  local  To  implement 
the  usual  leapfrog  FD-TD  time-stepping,  a  time  central-difference  equation  can  be  obtained  from 
(33)  by  treating  it  as  an  ordinary  differential  equation  in  time  with  the  right-hand  side  assumed  to 
be  a  constant  value  established  at  the  midpoint  of  the  time  step.  Note  that  this  procedure  keeps 
the  time  dependence  for  7(£).  As  shown  in  [3],  the  resulting  differential  equation  can  be 
accurately  time-stepped  for  a  variety  of  single  circuit  elements  (the  resistor,  capacitor,  inductor, 
and  diode). 

However,  for  the  case  of  one  of  the  ports  of  an  arbitrary  two-port  nonlinear  active  circuit 
connected  to  the  observation  point,  J  would  be  a  function  of  the  action  of  all  of  the  linear  and 
nonlinear  devices  comprising  the  circuit,  as  well  as  both  the  local  E  and  the  remote  E  sensed  at 
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the  other  port  of  the  circuit.  In  this  case,  the  simple  single-element  time-stepping  approach  of  [3] 
would  not  be  useful.  A  rationale  for  involving  the  general  circuit  simulator  SPICE  now  become 
apparent.  Because  SPICE  calculates  the  current  through  a  circuit  element  as  a  function  of  the 
voltage  across  the  device,  SPICE  could  provide  7(E)  for  a  circuit  device  simultaneously 
exposed  to  the  electromagnetic  field  and  embedded  in  a  general  electrical  network.  This  value  of 
J  could  be  used  in  the  left-hand  side  of  (33),  and  a  separate  time  central-difference  procedure 
could  be  done  to  provide  the  time-stepping  relation  for  E. 

However,  an  even  simpler  and  more  robust  solution  [30]  can  be  obtained  by  rewriting  (33)  as 


dV, 


-'grid  cell ' 


device  I  /  fV  .  ^  =  7  , 

'T  •■device V  'device /  ^ total 


dt 


(34) 


where 


device  the  voltage  across  (he  circuit  device,  and  therefore,  the  grid  cell; 

Cgridceii  =  eA/ A  is  the  FD-TD  grid-cell  capacitance,  where  A  is 

the  cross-section  area  of  the  grid  cell  and  A  is  its  height; 

4evice(^device)  =  A  7(E)  is  the  cunent  flowing  through  the  device, 

which  can  have  a  nonlinear  dependence  upon  ;  and 

is  the  known  total  current  flowing  through  the  grid  cell,  which  equals 
the  line  integral  of  the  H-field  vector  components  looping 
around  the  cell.  These  H  components  were  just  obtained  in  the 
magnetic  field  time-stepping  part  of  the  FD-TD  algorithm. 

Equation  (34)  can  be  represented  as  an  equivalent  circuit  consisting  of  a  current  source  in  parallel 
with  a  capacitor.  Thus,  instead  of  using  SPICE  just  to  determine  7(E),  SPICE  can  be  used  to 
time-step  (34)  directly  to  update  Vdevke*  which  then  yields  the  desired  updated  cell  electric  field 
by  the  simple  relation  E  =  Vjevice/A.  In  this  way,  the  lumped  element  can  be  an  arbitrarily  large 
SPICE  network  whose  description  can  be  contained  in  a  standard  SPICE  file.  Furthermore,  all  of 
the  extensive  device  models  in  SPICE  can  be  used  directly  in  the  FD-TD  simulation  without  the 
need  to  duplicate  the  model  development;  and  the  efficient  circuit  integration  methods  used  in 
SPICE  are  also  directly  available  without  user-implemented  integration  schemes. 

5.3  Reported  Applications 

This  approach  was  successfully  implemented  in  [30].  Initial  comparisons  were  made  between 
FD-TD/SPICE  and  pure  SPICE  simulations  for  one-port  circuits  placed  at  the  end  of  a  raicrosttip 
line.  These  included  individual  resistors,  capacitors,  inductors,  and  diodes,  and  then  progressed 
to  more  complicated  multi-element  circuits  consisting  of  several  series  and  parallel  components. 
All  FD-TD  data  were  in  excellent  agreement  with  the  appropriate  pure  SPICE  model. 

Reference  [30]  then  reported  the  modeling  of  a  more  interesting  nonlinear  two-port  network, 
a  VHF  tuned  amplifier  mounted  on  a  circuit  board.  The  overall  circuit  consisted  of  a  step- 
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function  voltage  source,  a  source  microstrip  line,  an  NPN  bipolar  junction  transistor  embedded 
with  LC  networks  for  base  and  collector  impedance  matching,  a  load  microstrip  line,  and  a 
resistive  load.  Excellent  agreement  on  the  order  of  1%  was  obtained  between  the  FD-TD/SPICE 
results  for  the  complicated  transistor  base-voltage  time  waveform  and  a  pure  SPICE  model  that 
was  constructed  to  account  for  the  physics  of  the  two  transmission  lines.  Similar  agreement  was 
observed  for  the  collector  voltage  waveform.  This  high  level  of  agreement  indicated  Aat  the 
FD-TD/SPICE  model  permitted  a  self-consistent  simulation  of  the  flow  of  electromagnetic  wave 
energy  in  both  directions  through  the  nonlinear  two-port  amplifier  network. 

Reference  [31]  reported  the  use  of  FD-TD/SPICE  to  perform  a  small-signal  analysis  of  a 
5.7-GHz  center-frequency  microwave  amplifier  consisting  of  a  three-terminal  packaged  FET,  its 
biasing  circuit,  the  associated  resistors  and  capacitors,  and  the  strip  line  interconnects.  Using  a 
modulated  Gaussian  pulse  excitation,  the  S-parameters  were  obtained  over  a  wide  frequency 
range  using  the  discrete  Fourier  transform  of  the  calculated  time  responses,  as  discussed  earlier. 
Solid  agreement  was  obtained  between  the  calculated  and  measured  5-parameters  over  the  range 
dc  to  9  GHz.  In  further  work  by  this  group,  a  power  amplifier  was  analyzed  for  large-signal 
excitation.  The  results  were  in  good  agreement  with  those  of  the  harmonic  balance  analysis  [32]. 

5.4  Potential  Applications 

The  hybrid  FD-TD/SPICE  tool  will  be  optimally  applied  when  the  speed  of  a  circuit  is  so  high 
and  its  physical  embedding  is  so  complex  that  modeling  electro-magnetic  wave  "artifacts"  is 
crucial.  Both  digital  and  analog  applications  are  expected: 

1.  Digital.  An  increasingly  wide  range  of  digital  applications  is  expected  as  clock  speeds 
approach  microwave  frequencies.  Present  engineering  practice  of  using  electromagnetic 
field  solvers  only  to  obtain  lumped-circuit  equivalences  for  subsequent  insertion  into 
conventional  SPICE  models  may  one  day  yield  to  direct  time-domain  modeling  of  logic 
operations  from  FD-TD/SPICE. 

2.  Analog.  Here  applications  will  include  analysis  of  linearity,  intermodulation,  harmonic 
generation,  and  conversion  efficiency  of  MMICs  embedded  in  compact  and  complex 
structures.  Another  category  of  applications  will  include  radiation,  especially  by  arrays 
of  surface  patch  antennas  excited  by  Gunn  diodes,  tunnel  diodes,  transistor  oscillators, 
and  MMICs  located  directly  at  the  antenna  [33]. 

FD-TD/SPICE  should  also  have  excellent  applicability  to  modeling  digital  or  analog  circuit  upset 
due  to  external  natural  and  manmade  electromagnetic  insults  such  as  lightning,  electromagnetic 
pulse,  and  high-power  microwaves. 
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MISSION 

OF 

AFRL/INFORMATION DIRECTORATE  (IF) 


The  advancement  and  application  of  information  systems  science  and 
technology  for  aerospace  command  and  control  and  its  transition  to  air, 
space,  and  ground  systems  to  meet  customer  needs  in  the  areas  of  Global 
Awareness,  Dynamic  Planning  and  Execution,  and  Global  Information 
Exchange  is  the  focus  of  this  AFRL  organization.  The  directorate’s  areas 
of  investigation  include  a  broad  spectrum  of  information  and  fusion, 
communication,  collaborative  environment  and  modeling  and  simulation, 
defensive  information  warfare,  and  intelligent  information  systems 


technologies. 


